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We investigate the evolution of multi-scale mechanical properties towards the macroscopic me- 
chanical instability in frictional granular media under compressive loading. Spatial correlations of 
shear stress redistribution following nucleating contact sliding events and shear strain localization 
are investigated. We report growing correlation lengths associated to both shear stress and shear 
strain fields that diverge simultaneously as approaching the transition to a dense flow regime. This 
shows that the transition from quasi static to dense flow regime can be interpreted as a critical 
phase transition. Our analysis does not show any characteristic shear band thickness potentially 
formed at the onset of instability. 
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The mechanical behavior of granular materials is of 53 
wide concern, from natural hazard in geological context 54 
to engineering application. However, the evolution of 55 
properties towards the flowing instability is still par- 56 
tially understood. Within packings of spheres, hetero- 57 
geneous distributions of contact forces on a scale much 58 
larger than the typical particle size pQ may control the 59 
mechanical response of the granular assembly, character- m 
ized by long range correlations and collective motions [2] . ei 
The concept of jamming [3] provides a powerful frame- 62 
work to analyse the onset of granular flows. Large ef- 63 
forts have been done in this direction by considering as- 64 
semblies of non frictional particles, which exibit jammed 65 
states resisting small stresses without irreversible defor- 66 
mation, whereas unjammed systems flow under any ap- 67 
plied stresses [3] . As granular materials play a fundamen- 68 
tal role in geophysical instabilities (e.g. granular gouges 69 
within fault zones, landslides), the concept of jamming 70 
might be a powerful tool to investigate these problems. 71 
However, two more ingredients are of primary importance 72 
in those cases. First, frictional granular materials have 73 
to be considered, leading to changes in the physics of 74 
jamming: the packing fraction is not a pertinent param- 75 
eter anymore and the transition from an unjammed to a 76 
jammed state is rather controlled by the non-rattler frac- 77 
tion [H [5] . Secondly, the loading conditions drastically ?s 
differ from the one commonly used to study jamming, 
since, rather than performed at constant volume, a con- 79 
fining pressure has to be considered. In that case, the so 
non-rattler fraction is not allowed to evolve freely and a ai 
percolating strong force newtork remains in the flowing 82 
phase, called the dense flow regime [6]. 83 



Here we consider compression tests under multi-axial 
loading: an increase of the axial stress a± is pre- 
scribed whereas the 2D confining pressure 03 is kept con- 
stant. These multi-axial conditions correspond to the 
generic case encountered in the Earth and the stress con- 
trol avoids stress relaxations and associated feedbacks, 
i.e. the stress and strain localization structures develop 
freely. Experiments [7] and theoretical works [5] have 
been conducted in this configuration on either continuous 
or discrete materials. In continuous materials, the macro- 
scopic instability as been first tackled by the use of bifur- 
cation theory, which considers a transition from an ho- 
mogeneous to an heterogeneous deformation field mate- 
rialized by the creation of a perennial macroscopic shear 
band spanning the whole sample [H [9] . Within granular- 
materials, either in experiments or simulations, comput- 
ing the deformation field over a large macroscopic strain 
window, those perennial shear bands appear and seem 
to show characteristic sizes [TUTU]. However, this vision 
is counterbalanced by the heteregenous and long range 
correlated kinematics of quasi-static granular flow [2] : at 
which temporal and spatial scales and at which stage of 
the loading can we consider the granular assembly to de- 
form homogeneously? What are the relevant key features 
of the deformation and stress fields associated to the on- 
set of instability? 

We use the Molecular Dynamics discrete element 
method to perform our simulations [llj . Two- 
dimensional granular assemblies made of frictional cir- 
cular grains are considered. The grains surfaces are 
uniformly distributed, setting the largest grain diame- 
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FIG. 1. Macroscopic behavior of a sample made of 10000 117 
grains during compressional testing. Top : Deviatoric stress iig 
r/<73 = (oi —crz)/cT3 (blue line) and macroscopic volume vari- 
ation AV/Vo (green line) as a function of axial deformation. 119 
Color dots correspond to locations where coarse-graining anal- 120 
yses of Figures [3|Top) and [4|Top) are performed. Bottom 121 
: Inertial number I as a function of the imposed deviatorici22 
stress t. I23 



conditions are used 1111 . Grains interact via linear elastic 

— „ 126 

laws and Coulomb friction when they are in contact |12j . 

Dense and highly coordinated initial packings, character- 
ized by a porosity r)i ~ 0.15 and a coordination number i2g 
Zi R3 4, are obtained by isotropic frictionless compression^ 
of dilute grains sets. We focus in this study on com- i 
pression tests, setting the particle friction to fi m i C ro — 1- 132 
The confining pressure 03 is sized by setting the contact^ 
stiffness k = k n /<r 3 equal to 1000 [11], where k n is the i34 
normal contact stiffness coefficient. The vertical stress <ji 

1 135 

is increased at constant rate. To characterize sample size 

136 

effects [TU] , we performed 320 simulations on 2500 grains 
samples, 80 simulations of 10000 grains and 20 simula- 
tions of 45000 grains. Most of the results presented here 138 
concern 10000 grains samples. Whatever the sample size, 139 
we used a stress increment Sa\ r equal to 1.10~ 6 x <j 3 at 140 



each discretisation time interval t r = 



X fe„ 



■/25, where 



mmin is the mass of the lighter grain. This leads to an™ 
inertial number / llj of the order of 10~ 6 — 10~ 5 at thei 4 4 
initial stage of the test (Figure [lj, which ensures quasi-145 
static conditions. lt6 
Figure [lja) shows a typical macroscopic stress-straini47 
curve. An initial contracting phase is observed, only be-«8 
cause of an elastic contact between particles [T3], untile 
a peak of contraction is reached, after which the sample^ 
dilates continuously. The deviatoric stress r = o\ — 173151 
reaches a plateau at values slightly larger than 2a 3 , ob-152 
tained for values of axial strain e\ of about 0.004. Be-153 
fore this specific point, while undergoing brutal increase^ 
during plastic avalanches, the inertial number / remainsiss 
lower than 10~ 4 (Figure [ijb)), which is often consideredise 
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FIG. 2. Left: Delaunay (top) and modified Voronoi (bot- 
tom) tesselations for a polydisperse granular material. Right: 
Coarse graining analysis on a 2500 grains sample. 



as the upper bound for quasi-static conditions [TTl [T3] . 
This region corresponds to the quasi-static regime. On 
the reverse, at values of r greater than r c = 203, that cor- 
responds to a macroscopic friction /i macro ~ 0.5, a brutal 
increase of / of several orders of magnitude is observed, 
reaching values of the order of 10~ 2 , which corresponds 
to a transition toward a dense flow regime, where inertia 
comes into play. At this transition, stress concentrations 
resulting from cooperative effects are expected, trigger- 
ing preferential weak zones where flow is favoured. In the 
case of our samples made of circular grains, this leads to 
a softening of the whole granular assembly and thus to a 
macroscopic friction much smaller than (J> m icro- Accord- 
ing to this, studying the spatial structure of both stress 
and strain fields is a key point to understand the mecha- 
nisms that generate the macroscopic instability. We thus 
focus, in this study, on the "dynamics" of the granular 
assembly, i.e. on the response of the granular assembly 
to a small stress increment in terms of associated stress 
concentration and strain localization structures that form 
during mechanical loading. 

We first characterize the spatial extent of stress con- 
centrations within the granular sample by means of a 
coarse graining analysis [TSl US] : an averaged shear 
stress rate < f > [17] is computed at different stages 
of mechanical testing (cf color dots on Figure [l]) over a 
time window T and over a broad range of spatial scales 
L, from the micro-scale corresponding to the scale of 
the mesh element, to the macroscale corresponding to 
sample size (Figure [2]). For a given assembly of grains, 
the components of the stress tensor are computed as 
a ij = V E c ( rC ~ r9 )fg [IS]' where V is the surface 
associated to the grain assembly, that corresponds to the 
sum of surfaces associated to each grain including the 
amount of porosity calculated using the modified Voronoi 
tesselation (Figure |2j), fg is the contact force exerted on 
grain g at contact c, r c is the position of c, and r 9 is 
the position of the center of g. Considering two succes- 
sive configurations, the stress rate tensor is obtained by 
differentiating the respective stress tensor components in 
time. The time resolution used is T 
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FIG. 3. Multi-scale analysis performed on the shear stress 
rate field f . A selection of corresponding fields is shown on 
the right. Top: < f > versus L for decreasing values of A. 
Configuration locations on the stress-strain curve are shown 
on Figure [l] The timescale T is used to compute values of 

< f >. The inset displays data collapse with respect to A 
(equation [TJ . We find v T = 1.3, C = 0.5 and 5=1. Bottom: 

< f > versus L at A = 0.005 when increasing the timescale 
t from t = T to t = 16T. For graphical convenience, all 
computed values have been normalized by the ones computed 
at the micro-scale. 



where N g is the number of grains of the sample. This 
timescale corresponds to the travel time of elastic waves 
through the granular assembly [10 . Results are shown 
on Figure [3] (Top) , from the early stages of biaxial test- 179 
ing up to t c = 2<7 3 . At the early stages of macroscopic 180 
deformation, a decrease of < f > with L is observed at 181 
small scales while for L-values larger than a crossover 182 
scale I* a plateau is observed. This means that shear 183 
stress rate fields are heterogeneous for I << I*, and ho- 184 
mogeneous ior I >> I*. Hence, I* is the associated corre- 185 
lation length [IB]. As macroscopic deformation proceeds, 186 
I* grows until reaching the entire size of the system at 187 
t c where a power law scaling < f >~ L~ Pt is observed, 188 
with p T = 0.38. The cut-off remaining on the scaling 189 
at A — > 0, where A = is the control parameter, 190 

is a finite size effect [TO]. These results suggest a pro- 191 
gressive structuring of the stress field as approaching the 192 
transition to the dense flow regime, associated to the di- 193 
vergence of the correlation length I*. It can be verified 194 
from a collapse analysis (inset of Figure ^ Top)) that I* 195 
diverges as L ' 
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FIG. 4. Multi-scale analysis performed on the incremental 
shear strain field Sy. A selection of corresponding fields is 
shown on the right. Top: < 5y > versus L for decreas- 
ing values of A towards the critical point. The deformation 
scale <5ei = 8e p = 1.10 -5 , consistant with the susceptibility 
analysis presented on Figure [5] is used to compute values of 
< <$7 >. The inset displays data collapse with respect to A 
(equation [l]. We find u-y = 1.3, C = 0.5 and 6 = 1. Bot- 
tom: < 8y > versus L at A = 0.005 when increasing <5ei from 
Sei = Se p to Sei = 32Se p . The inset displays data collapse 
with respect to £ = ^li_il£ "\y e g nc i an/ — q 4 p or graph- 
ical convenience, all computed values have been normalized 
by the ones computed at the micro-scale. 



where v T = 1.3 ±0.1 is the divergence exponent. Param- 
eters 8 and C characterize the finite size effect. A simi- 
lar analysis performed on other moments < r q > of the 
shear rate [TO] confirms the divergence of I* at A — > 0, 
and reveals their multi-fractality at the critical point. As 
stated previously, these particular features of the shear 
stress field are observed at a specific timescale t = T cor- 
responding to the clastic waves travel time. This multi- 
scale behaviour is no longer observed at larger timescales 
(Figure [^Bottom)), as a clear departure from power law 
is observed for t > T. Hence, the multi-scale properties 
of the shear stress rate field are only observed at the time 
of propagation of the elastic information throughout the 
sample. Beyond this time, a decorrelation is observed, 
explained by the superposition in time of several uncor- 
related events, consistant with a correlated stress struc- 
ture associated with little memory, limited to the elastic 
wave travel time. 

To study whether similar observations can be reported 
on the shear strain field, we consider a delaunay trian- 
gulation performed on the grain centers (Figure^]), af- 
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ter having removed the rattlers grains from the grain set, 
that allows to compute the partial derivatives at the mesh 
scale as e^- = l/2(dui/dxj+duj/dxi), where (1*1,^2) and 
(xx,X2) are respectively the incremental displacements 
and spatial coordinates of grain centers. Our coarse 
graining analysis is then performed by averaging par- 
tial derivatives at various scales, similarly to (author?) 
[HI [16]. An average shear strain rate [17] is thus obtained 
as a function of L. While not shown here, if one uses the 
constant timescale T to compute incremental displace- 
ments, a divergence of the correlation length associated 
with the shear strain rate field is not observed. Thus, 
the structure of the total strain field does not form si- 
multeanously to the stress field, by means of long-ranged 
elastic interactions. Intuitively, this would be the case 
if one would consider only the elastic component of the 



strain. Here, for N„ 



10000, 



10 stress increments 



are prescribed during the propagation time of an elastic 
wave throughout the sample. Hence, a multitude of con- 
tacts, in our case about 5% of the whole contact network, 
are then sliding, although elastic interactions didn't have 256 
time to travel across the entire sample, i.e. we are not 257 
in conditions of intinitly slow driving that could be ob- 
tained using event-driven algorithms [11] , Despite this, 
a progressive structuring of the shear strain field is ob-; 
served by considering constant macroscopic deformation; 
windows Sei = Se p = 1.10~ 5 to compute the scaling, 
of < <57 > (Figure [4]) . A divergence of the correlation, 
length I* similar to the one observed on the shear stress, 
rate field is obtained as approaching the transition to. 



the dense flow regime, as we find v 1 



1.3 from a„ 



231 collapse analysis (equation [T]) . When considering larger 267 

232 macroscopic deformation windows, the mutli-scale prop- 268 

233 erties of the deformation field are no longer observed, 269 

234 pointing out that Sei = Se p is a characteristic value. As 270 

235 the material softens when r increases, the corresponding 271 

236 time of integration decreases as the critical point is ap- 272 

237 proached. This corresponding time can either be smaller 273 

238 or larger than the elastic wave traveling time T, depend- 274 

239 ing on the imposed loading rate Sa\ r . However, whatever 275 

240 the loading rate considered and the size of the system, 276 

241 Sei — Se p remains equal to 1.10~ 5 , pointing out that a 277 

242 given amount of plastic activity has to operate in order 278 

243 to observe mutli-scale properties within the increment al 279 

244 shear strain field. As the correlation lengths I* and /* di- 280 

245 verge the same way, this plastic activity is probably the 281 

246 direct consequence of the progressive structuring of the 282 

247 stress field. 283 

248 An understanding of the characteristic value 5e p can 284 

249 be obtained from a four-point dynamic susceptibility 285 

250 X4 [HI [20] analysis on the inter-particles contact net- 28 6 

251 work. From a contact configuration that we refer as 287 

252 "initial" , selected at a value of axial deformation de- 288 

253 noted e™*, we compute the self-overlap order parame- 28 9 

254 ter Q e} „« (Sei) = ^ J2?=i ' 
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FIG. 5. Susceptibility analysis performed in the quasi-static 
phase on the sliding contacts taking part of the major (red 
curves) and minor (black curves) network. Top: Averaged 
self-overlap order parameter < Q(5ei) >. Bottom: Corre- 
sponding four-point average susceptibility < xa(S^i) >■ The 
vertical dashed line indicates the deformation value 5e p = 
1.10 , where a change of behavior is observed. 
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where N r is the number of29o 



255 contacts that are not sliding in the initial configuration 



and Wi is a step-function cutoff that equals 1 if no sliding 
event has been recorded on contact i over the whole de- 
formation window e™' — > e™* + Sei, and otherwise. 
259 The first two moments Q{5ei) =< Q £ «««(5ei) > and 
X i{5ei) = N c [< Q tT u{5eif > - < Q e *»«(*ei) > 2 ] 
of Q e imt(8ei) (calculated from sample-to-sample fluc- 
tuations) are then computed in the quasi-static region 
(t < r c ). Like this, we evaluate from an initial config- 
uration the number and the associated spatial hetero- 
geneity of sliding events nucleation as axial deformation 
increases. Figure [5] shows < Q(Sei) > and < x±($ e i) >j 
where < . > here means an average over all the val- 
ues of A (since no significant variation of Q(Sei) and 
X<i(8ei) is observed in the quasi-static region) computed 
by considering separately the major and minor force net- 
works. The major force network is defined by select- 
ing contact force greater than the average. By construc- 
tion, < Q(Sei) > is initialy equal to 1. As Sei increases, 

< Q(Sei) > decreases but never reaches 0, meaning that a 
considerable amount of contacts never slides. When con- 
sidering the whole test, about 35% (respectively 55%) 
of the contacts of the minor (respectively major) net- 
work didn't slide at the end of the test, meaning that the 
permanent deformation is extremely localized and that 
rigid bodies remain throughout the whole test [3T]. The 
value of < Xi(^ e l) > indicates, with respect to increas- 
ing axial deformation, the variability in the nucleation of 
new contact slidings. At low values of Sei, it increases 
as < X4(^ e i) >^ Se^ with (3 = 1.8, meaning that spa- 
tially correlated sites of contact sliding events are nucle- 
ating. As Sei exceeds the treshold value Se p = 1.10 -5 , 

< X4(^ £ i) > saturates, as all the spatially correlated con- 
tacts located close to the coulomb criteria, i.e. suscepti- 
ble to slide, have been destabilized. At this stage, only 
4% (respectively 3%) of contacts have slided at least one 
time in the minor (respectively major) network. 
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To conclude, incremental stress and strain fields are3« 
both characterized by a growing correlation length that 3 « 
diverges at the onset of macroscopic instability, which 342 
can therefore be identified as a critical point. A similar 343 
behavior has been reported in compressive failure of con- 
tinuous materials [TBI 122] • We interpret these stress and 
strain specific structures as resulting from dynamic stress 347 
redistributions induced by local dissipation of elastic en- 348 
ergy materialized by contact slidings. At macroscopic349 
instability, a local contact sliding event induces corre-350 
lated elastic stress perturbations up to the scale of the 351 
whole granular assembly. These features can only be ob- 352 
served when carefully examining characteristic timescales^ 
for stresses, and macroscopic deformation for strains con- 355 
centration processes. Thus, a long time memory, proba- 356 
bly associated to a slow evolution of the topology, existS357 
in the system, in relation with the growth of I* and Z* 358 
during mechanical loading towards the macroscopic in- 359 
stability, while the spatial multi-scales structure of stress 360 
and strain fields have respectively a short time and de- 
formation life. 363 

The last question that arises is to whether a limit in 364 
decreasing correlation length I* on the shear strain field 365 
is reached at A — > for values of Sei much greater than 
Se p , which would characterize the thickness of a peren- 368 
nial macroscopic shear band potentially formed at the 369 
onset of instability. To investigate this, we hypothesize37o 
that, close to the critical point (A — > 0), I* varies as 371 



where £ = 



<5ei — 5e 



L s is the square roo 



7 L*£ q t+C t S Se 

of the sample area and ct 7 is an exponent. This hypothe- 374 

sis is checked from a collapse analysis (inset of Figure [4]). 375 

We find a 1 — 0.7. This shows that I* keeps decreasing as 376 

the considered deformation window size Sei is increased, 377 

showing that the correlation length only depends on the 378 

value of <5ei and that no intrinsic scale of saturation, po- 

. 380 

tentially associated to a shear band thickness, can be 381 
identified at the onset of macroscopic instability. 382 
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